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Abstract: A scheme for rapidly and accurately computing solutions to boundary integral 
equations (BIEs) on rotationally symmetric surfaces in M 3 is presented. The scheme uses the 
Fourier transform to reduce the original BIE defined on a surface to a sequence of BIEs defined 
on a generating curve for the surface. It can handle loads that are not necessarily rotationally 
symmetric. Nystrom discretization is used to discretize the BIEs on the generating curve. The 
quadrature is a high-order Gaussian rule that is modified near the diagonal to retain high- 
order accuracy for singular kernels. The reduction in dimensionality, along with the use of 
high-order accurate quadratures, leads to small linear systems that can be inverted directly 
via, e.g., Gaussian elimination. This makes the scheme particularly fast in environments 
involving multiple right hand sides. It is demonstrated that for BIEs associated with the 
Laplace and Hclmholtz equations, the kernel in the reduced equations can be evaluated very 
rapidly by exploiting recursion relations for Legendre functions. Numerical examples illustrate 
the performance of the scheme; in particular, it is demonstrated that for a BIE associated 
with Laplace's equation on a surface discretized using 320 800 points, the set-up phase of the 
algorithm takes 1 minute on a standard laptop, and then solves can be executed in 0.5 seconds. 

1. Introduction 

The premise of the paper is that it is much easier to solve a boundary integral equation (BIE) 
defined on a curve in M 2 than one denned on a surface in IR 3 . With the development of high 
order accurate Nystrom discretization techniques [HI [2l [151 El E] > it has become possible to 
attain close to double precision accuracy in 2D using only a very moderate number of degrees 
of freedom. This opens up the possibility of solving a BIE on a rotationally symmetric surface 
with the same efficiency since such an equation can in principle be written as a sequence of BIEs 
defined on a generating curve. However, there is a practical obstacle: The kernels in the BIEs on 
the generating curve are given via Fourier integrals that cannot be evaluated analytically. The 
principal contribution of the present paper is to describe a set of fast methods for constructing 
approximations to these kernels. 

1.1. Problem formulation. This paper presents a numerical technique for solving boundary 
integral equations (BIEs) defined on axisymmetric surfaces in IR 3 . Specifically, we consider 
second kind Fredholm equations of the form 

(1) a{x) + j k(x,x')a(x')dA(x') = f(x), xeT, 

under two assumptions: First, that T is a surface in IR 3 obtained by rotating a curve 7 about 
an axis. Second, that the kernel k is invariant under rotation about the symmetry axis in the 
sense that 

(2) k(x,x') = k(9 - 8',r,z,r',z'), 

where (r, z, 6) and (r', z', 6') are cylindrical coordinates for x and x' , respectively, 

(3) x = (r cos6*, r sin 9, z), 

(4) x' = (r cos 6', r s'vuQ , z'), 

see Figure[TJ Under these assumptions, the equation 0, which is defined on the two-dimensional 
surface T, can via a Fourier transform in the azimuthal variable be recast as a sequence of 
equations defined on the one-dimensional curve 7. To be precise, letting a n , f n , and k n denote 
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the Fourier coefficients of a, /, and k, respectively (so that ([9]), (10), and (11) hold), the equation 
([!]) is equivalent to the sequence of equations 



(5) a n (r,z)+V2Tr k n (r,z,r',z')a n (r',z')r'dl(r',z') = f n (r,z), (r,z)ej, n € Z. 
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Whenever / can be represented with a moderate number of Fourier modes, the formula ^ 
provides an efficient technique for computing the corresponding modes of a. The conversion of 
@ to (|]) appears in, e.g., |20j, and is described in detail in Section [2j Note that the conversion 
procedure does not require the data function / to be rotationally symmetric. 

1.2. Applications and prior work. Equations like ([!]) arise in many areas of mathematical 
physics and engineering, commonly as reformulations of elliptic partial differential equations. 
Advantages of a BIE approach include a reduction in dimensionality, often a radical improve- 
ment in the conditioning of the mathematical equation to be solved, a natural way of handling 
problems defined on exterior domains, and a relative ease in implementing high-order discretiza- 
tion schemes, see, e.g., [3]. 

The observation that BIEs on rotationally symmetric surfaces can conveniently be solved by 
recasting them as a sequence of BIEs on a generating curve has previously been exploited in the 
context of stress analysis HJ, scattering dU E21 E3 123 , and potential theory [121 [TU1 1201 I2T] . 
Most of these approaches have relied on collocation or Galerkin discretizations and have generally 
used low-order accurate discretizations. 

1.3. Kernel evaluations. A complication of the axisymmetric formulation is the need to de- 
termine the kernels k n in ([5]). Each kernel k n is defined as a Fourier integral of the original 
kernel function k in the azimuthal variable 9 in Q, cf. ([8]), that cannot be evaluated analyt- 
ically, and would be too expensive to approximate via standard quadrature techniques. The 
FFT can be used to a accelerate the computation in certain regimes. For many points, however, 
the function which is to be transformed is sharply peaked, and the FFT would at these points 
yield inaccurate results. 

1.4. Principal contributions of present work. This paper resolves the difficulty of comput- 



ing the kernel functions k n described in Section 1.3 For the case of kernels associated with the 



Laplace equation, it provides analytic recursion relations that are valid precisely in the regions 
where the FFT loses accuracy. The kernels associated with the Helmholtz equation can then be 
obtained via a perturbation technique. 

The paper also describes a high-order Nystrom discretization of the BIEs ^ that provides 
far higher accuracy and speed than previously published methods. The discretization scheme 
converges fast enough that for simple generating curves, a relative accuracy of 10 -10 is obtained 
using as few as a hundred points, cf. Section [8| The rapid convergence of the discretization leads 
to linear systems of small size that can be solved directly via, e.g., Gaussian elimination, making 
the algorithm particularly effective in environments involving multiple right hand sides or when 
the linear system is challenging for iterative solvers (as happens for many scattering problems). 

Finally, the efficient techniques for evaluating the fundamental solutions to the Laplace and 
Helmholtz equations in an axisymmetric environment have applications beyond solving boundary 
integral equations, for details see Section [7j 

1.5. Asymptotic costs. To describe the asymptotic complexity of the method, we let Atot 
denote the total number of discretization points, and assume that as N tot grows, the number of 
Fourier modes required to resolve the solution scales proportionate to the number of discretiza- 
tion points required along the generating curve 7. Then the asymptotic cost of solving ([!]) for a 
single right-hand side / is 0(N^ ot ). If additional right hand sides are given, any subsequent solve 
requires only O(N^) operations. Numerical experiments presented in Section indicate that 
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the constants of proportionality in the two estimates are moderate. For instance, in a simulation 
with iVtot = 320 800, the scheme requires 1 minute for the first solve, and 0.49 seconds for each 
additional right hand side (on a standard laptop). 

Observe that since a high-order discretization scheme is used, even complicated geometries 
can be resolved to high accuracy with a moderate number A^ ot of points. 




Figure 1. The axisymmetric domain T generated by the curve 7. 



1.6. Outline. Section [2] provides details on the conversion of a BIE on a rotationally symmetric 
surface to a sequence of BIEs on a generating curve. Section [3] describes a high order methods 
for discretizing a BIE on a curve. Section [4] summarizes the algorithm and estimates its com- 
putational costs. Section [5] describes how to rapidly evaluate the kernels associated with the 
Laplace equation, and then Section [6] deals with the Helmholtz case. Section [7] describes other 
applications of the kernel evaluation techniques. Section [8] illustrates the performance of the 
proposed method via numerical experiments. 

2. Fourier Representation of BIE 

Consider the BIE under the assumptions on rotational symmetry stated in Section |1.1| 
(i.e. r is a rotationally symmetric surface generated by a curve 7 and that k is a rotationally 
symmetric kernel). Cylindrical coordinates (r,z,9) are introduced as specified in ([3]). We write 
r = 7 x T where T is the one-dimensional torus, usually parameterized by 6 £ (—it, it]. 



2.1, Separation of Variables. We define for n G Z the functions / n , cr n , and k n via 



(6) 
(7) 
(8) 



fn(r,z) 
cr n (r,z) 
k n (r,z,r',z') 



T \/2~7T 



-ind 



a(O t r,z)d0, 



-ind 



T v27T 



k(6,r,z,r',z')d9. 



Formulas Q, ([7]), and Q define / n , cr n , and fc n as the coefficients in the Fourier series of the 
functions /, a, and k about the azimuthal variable, 



(9) 

(10) 

(11) 



And 



f{x) = Y j ^=fn(r,z), 



e in(8-e>) 

k(x, x') = k{6 — 9' , r, z, r , z) = >^ / _ k n (r, z, r , z). 



2tt 



To determine the Fourier representation of (|lj) , we multiply the equation by e m ® /\/2tt and 
integrate 6 over T. Equation ^ can then be said to be equivalent to the sequence of equations 



(12) a n (r,z) + 



7XT 



-ind 



T v27T 



k(x,x')d6 



o(x')dA(x') = f n (r, z) 



n G Z. 



Invoking (11), we evaluate the bracketed factor in (12) as 

-,—inB r „—ind 



(13) 



T \/27r 



k{x,x')d0 



T \/27r 



k(e-e',r,z,r',z')de 



-ind' I e 



-in(e-B') 



T v 2tt 



k(9-e',r,z,r',z')de 



-in6' 



k n (r,z,r',z'). 



Inserting (13) into (12) and executing the integration of 0' over T, we find that ([!]) is equivalent 
to the sequence of equations 

(14) a n (r, z) + V2tt / k n (r,z,r',z')a n (r',z')r' dl(r',z') = f n (r,z), n G Z. 
For future reference, we define for n G Z the boundary integral operators /C n via 

(15) [K, n a n ](r,z) = V2tt k n (r,z,r',z')a n (r',z')r' dl(r',z'). 



Then equation (14) can be written 

(16) (/ + K n ) O n = f n , Tl G Z. 

When each operator / + K n is continuously invertible, we write the solution of ([!]) as 

ind 

(17) v(r, z,0) = J2 -7^1(1 + K n )- l fn}{r, z). 
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2.2. Truncation of the Fourier series. When evaluating the solution operator (17) in prac- 
tice, we will choose a truncation parameter N, and evaluate only the lowest 2N + 1 Fourier 
modes. If ./V is chosen so that the given function / is well-represented by its lowest 2N + 1 
Fourier modes, then in typical environments the solution obtained by truncating the sum (17) 
will also be accurate. To substantiate this claim, suppose that e is a given tolerance, and that 
N has been chosen so that 

N ind 

(18) II/- E ^rfn\\<z, 

at V27T 
n=—N 

We define an approximate solution via 

N e in6 

(19) O-approx = 2_] —={I + K, n )~f n . 

— ' V 27T 

n=—N 

From Parseval's identity, we then find that the error in the solution satisfies 

i 2 = j2 nu+K-r7nii 2 < j2 iK'+K-nm/nii 2 

\n\>N \n\>N 



< max \\(I + 1C n )- l \\ l > \\f n \\ 2 < max + K n y l \\ z e 
M>N J ^ \\n\>N 



It is typically the case that the kernel k(x,x') has enough smoothness that the Fourier modes 
k n (r, z,r' , z') decay as n — > oo. Then \\JC n \\ — > as n — > oo and + /C n ) _1 || —> 1. Thus, an 
accurate approximation of / leads to an approximation in a that is of the same order of accuracy. 
Figure [4] illustrates how fast this convergence is for Laplace's equation (note that in the case 
illustrated, the original equation is \a + Ka = f, and it is shown that | + /C n ) _1 | | — > 1/2). 

3. Nystrom discretization of BIEs on the generating curve 

We discretize the BIEs §5§ defined on the generating curve 7 using a Nystrom scheme. In 
describing the scheme, we keep the formulas uncluttered by discussing a generic integral equation 

o(x) + I k(x, x') cr(x') dl{x') = fix), x£j, 

J-y 

where 7 is a simple smooth curve in the plane, and A; is a weakly singular kernel function. 

3.1. Quadrature nodes. Consider a quadrature rule on 7 with nodes {xi}( =1 C 7 and weights 
{Wi}j =1 . In other words, for a sufficiently smooth function cp on 7, 



(20) ! ip(x)dl(x)^J2^ a 

J-y 



, (Xi)Wi. 

i=l 

For the experiments in this paper, we use a composite Gaussian rule with 10 points per panel. 
Such a rule admits for local refinement, and can easily be modified to accommodate contours 
with corners that are only piece-wise smooth. 



3.2. A simplistic Nystrom scheme. The Nystrom discretization of (20) corresponding to a 
quadrature with nodes {xi}\ =l takes the form 

1 

(21) Gi + ^aijUj = f(xi), i = 1, 2, 3, . . ., I, 
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where {ay}y = x are coefficients such that 

(22) I k(xi, x') a(x') dl(x') « ^ ay o-(s 3 -)> i = 1, 2, 3, 

The solution of (21) is a vector er = [ci]f =1 such that each erj is an approximation to cr(ccj). 

A simplistic way to construct coefficients ay so that (22) holds is to simply apply the rule 
(20) to each function x' i— > A;(ajj, a;') <r(a; / ) whence 

(23) ay = k(xi, Xj) Wj. 

This generally results in low accuracy since the kernel k(x, x') has a singularity at the diagonal. 
However, the formula (23) has the great advantage that constructing each ay costs no more 



than a kernel evaluation; we seek to preserve this property for as many elements as possible. 
3.3. High-order accurate Nystrom discretization. It is possible to construct a high-order 



(24) 
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discretization that preserves the simple formula (23) for the vast majority of coefficients ay 
|13j . The only coefficients that need to be modified are those for which the target point X{ is 
nearj^] the panel r holding Xj . In this conceptually constructed as follows: First 

map the pointwise values {&j} x er to their unique polynomial interpolant on r, then integrate 
this polynomial against the singular (or sharply peaked) functions x' \— > k(xi,x') using the 
quadratures of [15] Operationally, the end result is that ay is given by 

Yl™=i H x i> Vi,j, P ) v i,j,P if x % and x j are near ' 
k(xi, Xj) Wj if Xi and Xj are not near. 

In ( |24[ ), m is a small integer (roughly equal to the order of the Gaussian quadrature), the 
numbers Vy, p are coefficients that depend on 7 but not on k, and 2/y )P are points on 7 located 
in the same panel as Xj (in fact, j/y „ = yy p when j and j' belong to the same panel, so the 
number of kernel evaluations required is less than it seems). For details, see |13| . 

4. The full algorithm 

4.1. Overview. At this point, we have shown how to convert a BIE defined on an axisymmetric 
surface in M 3 to a sequence of equations defined on a curve in R 2 (Section [2]), and then how 
to discretize each of these reduced equations (Section [3]) . Putting the components together, we 
obtain the following algorithm for solving ([I]) to within some preset tolerance e: 

i) Given /, determine a truncation parameter N such that || / — Yln=-N ^73= /"II — £ - 

ii) Fix a quadrature rule for 7 with nodes {(n, £j)}f =1 C 7 and form for each Fourier mode 
n = —N, —N +1, — A^+2, . . . , N the corresponding Nystrom discretization as described 
in Section|3j The number of nodes I must be picked to meet the computational tolerance 
e. Denote the resulting coefficient matrices {A^ nS) }^ = _ N . 

iii) Evaluate via the FFT the terms {/ n (n, Zi)}% = _ N (as defined by Q) for i = 1, 2, 3, ... , I. 

iv) Solve the equation (I + A( n )) a n = f n for n = -N, -N + 1, —N + 2, . . . , N. 



v) Construct <7 a p prox using formula (19) evaluated via the FFT 



The construction of the matrices A^ n ^ in Step ii can be accelerated using the FFT (as described 



in Section 4.2), but even with such acceleration, it is typically by a wide margin the most 



expensive part of the algorithm. However, this step needs to be performed only once for any 



"*To be precise, we say that a point Xi is near a panel r if Xi is located inside a circle that is concentric to the 
smallest circle enclosing r, but of twice the radius. 
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given geometry. The method therefore becomes particularly efficient when ([!]) needs to be solved 
for a sequence of right-hand sides. In this case, it may be worth the cost to pre-compute the 
inverse (or LU-factorization) of each matrix I + A( n ). 

4.2. Cost of computing the coefficient matrices. For each of the 27V + 1 Fourier modes, we 
need to construct an I x I matrix A( n ) with entries af 1 ) . These entries are derived from the kernel 
functions k n defined by (|8j). Note that whenever (r, z) ^ (r', z'), the function 9 h-> k(9, r, z, r' , z') 
is C°°, but that as (r', z') — > (r, z) it develops a progressively sharper peak around 9 = 0. 



For two nodes (rj, Zi) and (rj, Zj) that are "not near" (in the sense defined in Section 3.3) the 
matrix entries are given by the formula 

(25) aVV = k n (ri, Zi,rj,Zj) wj 

where k n is given by ([8]). Using the FFT, all 2N + 1 entries can be evaluated at once in 
0(N log N) operations. The FFT implicitly evaluates the integrals Q via a trapezoidal rule 
which is highly accurate since the points (r*j, Zj) and (rj, Zj) are well-separated on the curve 7. 

For two nodes (r^z,) and (rj,Zj) that are not well-separated on 7, evaluating af^ is dicier. 



The first complication is that we must now use the corrected formula, cf. (24), 

m 

(26) «S? = E fc 



I Ti , Zi 1 



The second complication is that the FFT acceleration for computing the kernels {k n }^ = _ N 
jointly no longer works since the integrand in ^ is too peaked for the simplistic trapezoidal 
rule implicit in the FFT. Fortunately, it turns out that for the kernels we are most interested 
in (the single and double layer kernels associated with the Laplace and Helmholtz equations), 
the sequence {k n } n __ N can be evaluated very efficiently via certain recurrence relations as 
described in Sections [5] and [6] (for the Laplace and Helmholtz equations, respectively). Happily, 
the analytic formulas are stable precisely in the region where the FFT becomes inaccurate. 



4.3. Computational Costs. The asymptotic cost of the algorithm described in Section 4.1 
will be expressed in terms of the number N of Fourier modes required, and the number I of 
discretization points required along 7. The total cost can be split into three components: 

(1) Cost of forming the matrices {A^}^ = _ N : We need to form 2N + 1 matrices, each of 
size I x I. For 0(I 2 ) entries in each matrix, the formula (25 ) applies and using the FFT, 



all 2N + 1 entries can be computed at once at cost 0(N log N). For O(I) entries close 



to the diagonal, the formula (26) applies, and all the 2N + 1 entries can be computed 
at once at cost O(N) using the recursion relations in Sections [5] and [6| The total cost of 
this step is therefore 0(I 2 N log N). 

(2) Cost of transforming functions from physical space to Fourier space and back: The 
boundary data / must be decomposed into Fourier modes {fn}n=-N> an< ^ after the 
linear systems (I + A^)a n = f n have been solved, the Fourier modes {c n }^ = _ N must 
be transformed back to physical space. The asymptotic cost is 0(IN log(A^)). 

(3) Cost of solving the linear systems (I + fv- n ')a n = f n : Using a direct solver such as 
Gaussian elimination, the asymptotic cost is 0(I 3 N). 

We make some practical observations: 

• The cost of executing FFTs is minimal and is dwarfed by the remaining costs. 
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• The scheme is highly efficient in situations where the same equation needs to be solved 
for a sequence of different right hand sides. In this situation, one factors the matrices 
I + A( n ) once at cost 0(I 3 N), and then the cost of processing an additional right hand 
side is only 0(I 2 N + IN log N) with a very small constant of proportionality. 

• To elucidate the computational costs, let us express them in terms of the total number 
of discretization points A^ot under the simplifying assumption that / ~ N. Since A^ot = 
IN, we find I ~ N^ and N ~ N^. Then: 

Cost of setting up linear systems: 0(iV tot logiVtot) 
Cost of the first solve: 0{N% ot ) 

Cost of subsequent solves: 0(iV tot ) 

We observe that even though the asymptotic cost of forming the linear systems is less 
than the cost of factoring the matrices, the set-up phase tends to dominate unless Atot 
is large. Moreover, the 0(-/V tot ) cost of the subsequent solves has a very small constant 
of proportionality. 

• The high order discretization employed achieves high accuracy with a small number of 
points. In practical terms, this means that despite the 0(iV t 2 ot ) scaling, the scheme is 
very fast even for moderately complicated geometries. 

• The system matrices I + A^ n ^ often have internal structure that allow them to be inverted 
using "fast methods" such as, e.g., those in [18]. The cost of inversion and application 
can in fact be accelerated to near optimal complexity. 

5. Accelerations for the Single and Double Layer Kernels Associated with 

Laplace's Equation 

This section describes an efficient technique based on recursion relations for evaluating the 
kernel k n , cf. ([8]), when k is either the single or double layer kernel associated with Laplace's 
equation. 

5.1. The Double Layer Kernels of Laplace's Equation. Let D C M 3 be a bounded domain 
whose boundary is given by a smooth surface T, let E = D c denote the domain exterior to D, 
and let n be the outward unit normal to D. Consider the interior and exterior Dirichlet problems 
of potential theory [TT] , 

(27) An = in D, u = f on T, (interior Dirichlet problem) 

(28) An = in E, u = f on T. (exterior Dirichlet problem) 



The solutions to (27) and (|28|) can be written in the respective forms 



(29) u(x) = f n (fy( x f\ ^ x ') d A(x'), xeD, 

Jy 4tt\x — x'\ 6 

, \ i \ I f n(x') ■ (x — x') 1 \ . , . . _ 
30 u(x)= / 7^ + ^ :)a{x')dA{x'), x£E, 

where a is a boundary charge distribution that can be determined using the boundary conditions. 
The point xq can be placed at any suitable location in D. The resulting equations are 

(31) -\a(x) + / n( f]^ X ~f <r(x') dA(x>) = f(x), 

2 Jy vk\x — x'l^ 

, x 1 , n f f n(x') ■ (x — x') 1 \ n / \ 

32 -^<j(x) + / \ V + —. r a{x') dA(x') = f(x), 

2 J r V 4:it\x - x'\ 3 4tt\x-x \J 
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where x G T in (31) and (32). 

Remark 1. There are other integral formulations for the solution to Laplace's equation. The 
double layer formulation presented here is a good choice in that it provides an integral operator 
that leads to well conditioned linear systems. However, the methodology of this chapter is equally 
applicable to single-layer formulations that lead to first kind Fredholm BIEs. 

5.2. Separation of Variables. Using the procedure given in Section [2j if T = 7 x T, then 
(27) and (28) can be recast as a series of BIEs defined along 7. We express n in cylindrical 
coordinates as 

Further, 



n(x') = (n r > cos 9' , n r i sin 9 , n z 



x — x 



/|2 



r' cos9') 2 + (rsin( 



r'sin^') 2 + (z — z') 2 



(r cos ( 

r 2 + (V) 2 — 2rr'(sin 9 sin 9' + cos 9 cos 9') + (z — z') 2 
r 2 + (r') 2 - 2rr' cos{9 - 9') + (z - z') 2 



and 



n(x ) • (x — x') = (n r i cos 9' , n r i sin 9' , n z /) • (r cos 9 — r' cos 0', r sin0 — r' sin#', z — z) 

= ?vr(sin 9 sin 9' + cos 9 cos 0') — ra r T + n z i (z — z) 

= n r '(rcos(6 — 9') — r) + n z i{z — z). 

Then for a point x' S T, the kernel of the internal Dirichlet problem can be expanded as 

n(s') • {x - x') = 1 e in(e-0') d (t) (r 
4ir\x - x'\ 3 -/o^ ^ n 



1 i\ 

z,r ,z) 



no 



where 



V 327T 3 



?V (r cos 6* — r') + n z / (z — z') 
r 2 + (y)2 _ 2 rr / cos 6) + [z - z ') 2 ) 3 / 2 



d9. 



Similarly, the kernel of the external Dirichlet problem can be written as 



n(x') • (x — x') 



4ir\x — x 



'13 



+ 



1 



with 



4 e) (wV) 



1 



4ir\x — xq\ 



-ind 



\/327r 3 J t 



+ 



^^e^')d(f)(r,z,r',z'), 

?V (r cos 9 — r') + n z / (z — z') 
.2 + ( r /)2 _ 2rr ' cos # + ( z - z') 2 ) 3 / 2 

d9, 



+ 



(r 2 + Tq — 2rr$ cos + (z — zq) 2 ) 1 / 2 _ 

where 2:0 has been written in cylindrical coordinates as (ro cos(#o)> r o sm ($o)i z o)- With the 
expansions of the kernels available, the procedure described in Section [4] can be used to solve 
(31) and (32) by solving 

(33) 

and 
(34) 

J-y 

respectively for n = —N, —N + 1, ... ,N. Note that the kernels d$ and d$ contain a log- 



cr n (r, z) + V2tt / d$(r, r', z, z')a n (r', z') r' dl(r', z') = f n (r, z) 

a n (r,z) + V2tt [ d^ (r,r , z, z)a n {r' , z) r dl(r' , z) = f n (r, z), 

n = -N,-N 
singularity as (r' , z') — > (r, z). 
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(i) (e) 

5.3. Evaluation of Kernels. The values of d n and d n for n 



-N, -N + 1, . . . ,N need to 



be computed efficiently and with high accuracy to construct the Nystrom discretization of (|33|) 

egr, 

interval [— vr, it] . Therefore, d$ can be written as 



and ( 34 ) . Note that the integrands of d n and d n are real valued and even functions on the 



(35) d^(r,z,r',z') 



1 



V327T 3 Jt 



n r / (r cos t — r') + n z i (z 



(r 2 + (r') 2 - 2rr' cos t + (z - z') 2 ) 3 / 2 



cos(nt) dt. 



(e) 

Note that dn can be written in a similar form. 

This integrand is oscillatory and increasingly peaked at the origin as (r', z') approaches (r, z). 
As long as r' and r as well as z' and z are well separated, the integrand does not experience 
peaks near the origin, and as mentioned before, the FFT provides a fast and accurate way for 
calculating dn and dn ^ ■ 

In regimes where the integrand is peaked, the FFT no longer provides a means of evaluating 
d$ and dff with the desired accuracy. One possible solution to this issue is applying adaptive 
quadrature to fully resolve the peak. However, this must be done for each value of n required 
and becomes prohibitively expensive if N is large. 

Fortunately, an analytical solution to (35) exists. As noted in jS], the single-layer kernel can 
be expanded with respect to the azimuthal variable as 



s(x, x') 



1 



1 



Airlx — x' 



4ir{r 2 + (r') 2 - 2rr' cos((9 - 9') + (z - z') 2 ) 1 / 2 
27T ~L 



where 



s n (r,z,r',z') 



cos(nt) 



v / 32~^ J T (r 2 + (r') 2 - 2rr' cos(t) + (z - z 1 ) 2 ) 1 / 2 



dt 



1 



cos(nt) 



V8ir 3 rr' Jt \/8{x ~ cos(t)) 
=== Qn-l/2(x) 



dt 



Qn-i/2 is the half- integer degree Legendre function of the second kind, and 

„/\2 



X 



r 2 + (r') 2 + (z 
2rr' 



To find an analytical form for ( 35 ) , first note that in cylindrical coordinates the double-layer 
kernel can be written in terms of the single-layer kernel, 



n(x 



(x — x 



n r i(r cos(# 



r') + n z i(z 



Air\x — x 



r\3 



4vr(r 2 + (r') 2 - 2rr' cos(fl - 9') + (z - z') 2 ) 3 / 2 



1 

4tt 



+ 7W- 



d_ 
dr' 




( r 2 + ( r /)2 _ 2rr> cos(6> - 9') + (z - z') 2 ) 1 / 2 



+ 



dz> \ (r 2 + (r') 2 - 2rr> cos(6> - 9>) + (z - z') 2 ) 1 / 2 
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(i) 

The coefficients of the Fourier series expansion of the double-layer kernel are then given by dn , 
which can be written using the previous equation as 

d (o (r Z1 / j)- n , [ JL( cos ( nt ) ^ dt+ 

a n \r,z,r,z) n r ^ ^ y + (r , )2 _ ^ + ^ _ ^ )2)) i/2 J at+ 

i d ( cos(nt) \ 

" " ' ' dz 1 V (327r 3 (r 2 + (r') 2 - 2rr' cos(t) + (z - z') 2 )) 1 / 2 ) * 



-»<■'— ! -^==^2n-l/2(x)) + n -'^7 ( 7==^ gn-l/2(x) 



1 



V87r 3 rr' 



, 0Qn-i/ 2 (x) 0X Qn-i/ 2 (x)\ 9Qn-i/2(x)a X 

»V ^ 7T^ 7T-. + n, ■ 



dx dr' 2r' J dx dz 



To utilize this form of d$ , set /i = y^r[ an< ^ no * e ^ na ^ 

dx _ {r' f -r 2 -(z- z') 2 
dr' 2r(r') 2 
dx z' — z 
dz' rr' ' 
S-i/ 2 (x) = nK(ji\ 

Qi/ 2 (x) = - ^2(x + i)^(m), 

Q-n-l/ 2 (x) = Qn-1/ 2 (X), 

n — 1 2n — 3 

S n -1/2(X) = 4 — _ Xgn-3/2(Xj ~ j^Qn-5/2(Xj, 

dQ n -i/2(x) In- 1 , . 
~~ 2(x 2 - 1) ^ n - 1 / 2 ~ y ™-3/2j i 

where K and 22 are the complete elliptic integrals of the first and second kinds, respectively. The 
first two relations follow immediately from the definition of x an d the relations for the Legendre 
functions of the second kind can be found in pp. With these relations in hand, the calculation 
of d y n> for n = —N, —N + 1,...,N can be done accurately and efficiently when r and r as well 
as z and z are in close proximity. The calculation of d n can be done analogously. 

Remark 2. Note that the forward recursion relation for the Legendre functions Q n _i/ 2 (x) w 
unstable when x > 1- practice, the instability is mild when x is near 1 and the recursion 
relation can still be employed to accurately compute values in this regime. Additionally, if stabil- 
ity becomes an issue, Miller's algorithm [TU] can be used to calculate the values of the Legendre 
functions using the backwards recursion relation, which is stable for x > 1 • 

6. Fast Kernel Evaluation for the Helmholtz Equation 

Section [5] describes how to efficiently evaluate the kernels k n as defined by Q for kernels 
associated with Laplace's equation. This section generalizes these methods to a broad class of 
kernels that includes the single and double layer kernels associated with the Helmholtz equation. 

6.1. Rapid Kernel Calculation via Convolution. Consider a kernel of the form 

(36) f(x, x') = s(x, x') g(x, x'), 

where 

f '\ 1 

S{X,X ) = — 77 

4tt\x-x'\ 
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is the single layer kernel of Laplace's equation and g(x,x') is a smooth function for all x,x' G 
M 3 . Common examples of kernels that take this form include the fundamental solution of the 
Helmholtz equation and screened Coulomb (Yukawa) potentials. 
Letting 

x = (r cos#, r sin (9, z), 
x = (r cos 6', r sin (9', z), 



we are interested in calculating the Fourier expansion of ( |36| ) in terms of the azimuthal variable. 
When g(x,x') = 1 (the Laplace kernel), we know how to rapidly compute these Fourier coef- 
ficients rapidly and efficiently. However, when g takes a nontrivial form, this is not generally 



true; there is no known analytical formula for calculating the Fourier coefficients of (36). 



We will now describe an efficient technique for calculating the the Fourier coefficients of ( 36 ) , 
when the function g is sufficiently smooth. For a fixed value of (r, z) and (r', z'), the functions s 
and g are periodic in the azimuthal variable over the interval T. Dropping the dependence of s 
and g on (r, z) and (r', z') for notational clarity, we define t = 9 — 6' € T and the Fourier series 
expansions of s and g as 



(37) 

(38) 

where 
(39) 

(40) 



s(t) 



Ant 



n int 



-int 



T v27T 



s(t) dt, 



9n 



n —int 



T v27T 



g(t) dt. 



The values given by (39) can be calculated as described in Section [5j while the values given by 
(40) can be rapidly and accurately computed using the FFT. 



Assuming the Fourier series defined by (37) and (|38|) are uniformly convergent, we find 

1 



fn 



2tt Jj 



/ s(t)g(t)e mt dt = \2s k g n - k = [s k * g k ](n), 



where s k * g k is the discrete convolution of the sequences defined by ( 39 ) and ( 40 ) . 



In a practical setting, the Fourier series are truncated to finite length. Assuming that we 
have kept — N, —N + 1, . . . , N terms, directly calculating the convolution would require 0(N 2 ) 
operations. Fortunately, this computation can be accelerated to 0(N log N) operations by 
employing the discrete convolution theorem and the FFT [6] . 

Letting T> denote the discrete Fourier transform (DFT), the discrete convolution theorem 
states that the convolution of two periodic sequences {a n } and {b n } is related by 

D{a n * b n } k = aA k B k , 

where {A n } = V{a n }, {B n } = V{b n }, and a is a known constant depending upon the length 
of the periodic sequences and the definition of the DFT. Thus, we can rapidly calculate the 
convolution of two periodic sequences by taking the FFT of each sequence, computing the 
pointwise product of the result, and then applying the inverse FFT. 

Of course, the sequences that we need to convolute are not periodic. Applying the discrete 



convolution to the sequences defined by (39) and (40) will not be exact, but the error incurred 
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will be small assuming that the Fourier coefficients decay rapidly and that ./V is large enough. 
To see this, assume that 



N 



»(*) = E 



n int 



n=-N 
N 



27T 



Ant 



9{t)= ^ 7^' 

n=—N 



Then the exact Fourier representation of / can be found by taking the product of these two 
series, which will be of length AN +1. As is well known, the coefficients of this product is given by 
the discrete convolution of the sequences containing the coefficients of the two series, and these 
sequences must first be padded with 2N zeros. Thus, we can effectively calculate the Fourier 



coefficients of the function given by (36) by calculating 2N + 1 Fourier coefficients of s and g, 



padding these sequences with zeros, calculating the discrete convolution of these sequences, and 
truncating the resulting sequence. In practice, padding may not even be required if the Fourier 
coefficients of s and g decay sufficiently fast. 

Note that the procedure described in this section is quite general. The azimuthal Fourier 
coefficients of many kernels that can be represented as the product of a singular function and 
a smooth function can be found, assuming that there is an accurate technique for determining 
the coefficients of the singular function. 

6.2. Application to the Helmholtz Equation. In this section, we will apply the fast kernel 



calculation technique described in Section 6.1 to the exterior Dirichlet problem for the Helmholtz 
equation. Let D C M 3 be a bounded domain whose boundary is given by a smooth surface T, 
let E = D c denote the domain exterior to D, and let n and be the outward unit normal to D. 
The partial differential equation representing this problem is given by 

(41) An + k 2 u = in E, u = f on T, 

where k > is the wavenumber, and u satisfies the Sommerfeld radiation condition 



(42) 



lim r [ — — 

? — s>oo v or 



iku 



0. 



where r = \x\ and the limit holds uniformly in all directions x/\x\. Let the single and double 
layer potentials for the Helmholtz equation be given by 



(43) 
(44) 



cp[x,x , 

dcf)(x, x') 
dn(x') 



Air\x — x'\ 
n{x') • [x — x' 



47T \x 



X 
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(single layer) 



(l-ik\x-x'\) 



ik\x—x'\ 



(double layer) 



The solution to (41) can be written in terms of the double layer potential, 

r d(f>(x,x') 
l r dn(x') 



u{x) 



a(x')dA(x'), x £ E, 



where a is a boundary charge density that can be determined using the boundary conditions. 
The resulting boundary integral equation is given by 

I f d(f)(x,x') 



(45) 



a(x')dA(x') = f(x). 



As is well known, (45) is not always uniquely solvable, even though (41) is uniquely solvable for 



all k > 0. A common solution to this is to represent the solution to (41) as a combined single 
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and double layer potential, 

u(x) = J ~ i " <Kx, x')\ a(x') dA(x'), xeE, 

where v > 0. We have freedom in choosing v, see, e.g., [3 [TO], for some analysis on this choice. 
The boundary integral equation we need to solve is 



(46) 



d<j){x, x') 
dn(x') 



iv<t>(x,x') \ a(x')dA(x') = f(x). 



7. Fast evaluation of fundamental solutions in cylindrical coordinates 

The techniques for kernel evaluations described in Sections |4.2[ [5j and [6] are useful not only 
for solving BIEs, but for solving the Laplace and Helmholtz equations in a variety of contexts 
where cylindrical coordinates are effective. To illustrate, observe that the free space equation 

(47) - Au(x) - k 2 u(x) = f(x), x E R 3 
has the solution 

(48) u(x)= [ cj) {k Xx,x') f{x')dA{x'). 
where (f)^ is the fundamental solution 

#)( 5 



x e 



^ik\x\ 



X 



4tt\x\ 



In cylindrical coordinates, we write (48) as 

oo 

E 



(49) 



:v e in9 



U(X) 



27T J H 



^\r,z,r',z')f n (r',z')dA(r',z') 



where H = {(r,z) £ R 2 : r > 0} is a half-plane, and where u n , f n , and (f^f 1 are the Fourier 
coefficients defined by 



u(x) = 

n=—c 
oo 

f(*) = E 



v e ind 



u n (r,z), 



fn(r,z), 



oo in{8—8') 

<P^(x,x')= e -^^<Pl k) (r,z,r',z'). 



2tt 



The kernel (j)^ in (49) can be evaluated efficiently using the techniques of Sections 4.2[ 5] and 
ro] If / has a rapidly convergent Fourier series, and if "fast" summation (e.g. the Fast Multipole 



Method) is used to evaluate the integrals in (49), then very efficient solvers result. 



More generally, we observe that the equation (|47|) can be expressed 
(50) 



d 2 u n 1 du n d 2 u % 



d 2 7 



r dr d 2 z 



+ 



n 



k I Un — fni 



n E 



(k) , 

and that the function (fy n is the Green's function of (50). 
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J 



(a) 





(b) 





(c) 

Figure 2. Domains used in numerical examples. All items are rotated about 
the vertical axis, (a) An ellipse, (b) A wavy block, (c) A starfish torus. 



8. Numerical Results 

This section describes several numerical experiments performed to assess the efficiency and 
accuracy of the numerical scheme outlined in Section |4.1| The geometries investigated are 
described in Figure [2j The generating curves were parameterized by arc length, and split into 
iVp panels of equal length. A 10-point Gaussian quadrature has been used along each panel, 
with the modified quadratures of [15] used to handle the integrable singularities in the kernel. 
The algorithm was implemented in FORTRAN, using BLAS, LAPACK, and the FFT library 
provided by Intel's MKL library. All numerical experiments in this section have been carried 
out on a Macbook Pro with a 2.4 GHz Intel Core 2 Duo and 4GB of RAM. 



8.1. Laplace's equation. We solved Laplace's equation in the domain interior to the surfaces 
shown in Figure [2] The solution was represented via the double layer Ansatz ( 29 ) leading to 
the BI E (|3lD . The kernels d n defined via ( 35 ) were evaluated using the techniques described in 
Section 5.3 Tn this case 

A (n) = A (-n) > and 

so we need only to invert N + 1 matrices. Further, 
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Figure 3 . Timings of the algorithm as the number of degrees of freedom iV tot = 
I(2N + 1) increases. The timings reported here are for the case I ~ 2N + 1. The 
numbers in parentheses provide estimates of the asymptotic complexity, i.e. the 
best fit to a curve T = C . 



the FFT used here is complex-valued, and a real-valued FFT would yield a significant decrease 
in computation time. 

To investigate the speed of the proposed method, we solved a sequence of problems on the 
domain in Figure [2^a). The timing results are given in Table [T] The reported results include: 

iVp the number of panels used to discretize the contour (each panel has I/Np nodes) 

N the Fourier truncation parameter (we keep 2N + 1 modes) 

T mat time to construct the linear systems (utilizing the recursion relation) 

Ti nv time to invert the linear systems 

Tfft time to Fourier transform the right hand side and the solution 
T app i y time to apply the inverse to the right hand side 

The most expensive component of the calculation is the kernel evaluation required to form 
the coefficient matrices . Table [2] compares the use of the recursion relation in evaluating 
the kernel when it is near-singular to using an adaptive Gaussian quadrature. The efficiency of 
the recursion relation is evident. 

Figure [3] plots the time to construct the linear systems as the number of degrees of freedom 
Ntot = I(2N + 1) increases, for the case when / ~ 2N + 1. The estimated asymptotic costs 



given in this Figure match well with the estimates derived in Section 4.3. It is also clear that 



as iVtot grows, the cost of inversion will eventually dominate. We remark that the asymptotic 
scaling of this cost can be lowered by using fast techniques for the inversion of boundary integral 
operators, but that little gain would be achieved for the problem sizes considered here. 

We observe that the largest problem reported in Table [I] involves 320 800 degrees of freedom. 
The method requires 1 minute of pre-computation for this example, and is then capable of 
computing a solution u from a given data function / in 0.49 seconds. 
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N P 


2iV + 1 


^rriat 


^inv 


T st 


^apply 


5 


25 


1.70E-02 


1.42E-03 


7.81E-05 


3.83E-05 


10 


25 


3.64E-02 


6.15E-03 


1.68E-04 


2.66E-04 


20 


25 


9.73E-02 


3.52E-02 


3.69E-04 


2.10E-03 


40 


25 


3.09E-01 


2.35E-01 


6.69E-04 


4.82E-03 


80 


25 


1.20E+00 


1.88E+00 


1.36E-03 


2.85E-02 


5 


51 


2.83E-02 


3.02E-03 


2.38E-04 


1.13E-04 


10 


51 


6.71E-02 


1.23E-02 


4.48E-04 


6.58E-04 


20 


51 


2.02E-01 


7.42E-02 


9.17E-04 


2.63E-03 


40 


51 


7.28E-01 


4.94E-01 


1.92E-03 


1.03E-02 


80 


51 


3.07E+00 


3.73E+00 


3.59E-03 


6.17E-02 


5 


101 


5.08E-02 


5.48E-03 


7.27E-04 


2.38E-04 


10 


101 


1.35E-01 


2.27E-02 


1.41E-03 


1.33E-03 


20 


101 


4.39E-01 


1.32E-01 


2.73E-03 


4.60E-03 


40 


101 


1.98E+00 


1.04E+00 


6.20E-03 


2.14E-02 


80 


101 


7.13E+00 


7.04E+00 


1.12E-02 


1.12E-01 


5 


201 


1.07E-01 


1.06E-02 


1.80E-03 


6.36E-04 


10 


201 


3.33E-01 


4.96E-02 


3.83E-03 


2.79E-03 


20 


201 


1.12E+00 


2.73E-01 


7.05E-03 


9.46E-03 


40 


201 


4.63E+00 


1.88E+00 


1.46E-02 


4.15E-02 


80 


201 


1.71E+01 


1.41E+01 


2.93E-02 


2.15E-01 


5 


401 


1.87E-01 


2.15E-02 


3.43E-03 


1.48E-03 


10 


401 


5.85E-01 


9.51E-02 


6.59E-03 


5.05E-03 


20 


401 


2.15E+00 


5.42E-01 


1.33E-02 


1.91E-02 


40 


401 


8.40E+00 


3.71E+00 


2.78E-02 


7.98E-02 


80 


401 


3.15E+01 


2.83E+01 


5.56E-02 


4.34E-01 



Table 1. Timing results in seconds performed for the domain given in Figure 
2[a) for the interior Dirichlet problem. 



2N + 1 


Composite Quadrature 


Recursion Relation 


25 


1.9 


0.017 


50 


3.1 


0.028 


100 


6.6 


0.051 


200 


18.9 


0.107 



Table 2. Timing comparison in seconds for constructing the matrices (I + A^) 
using composite Gaussian quadrature and the recursion relation described in 
Section 5.3 to evaluate k n for diagonal and near diagonal blocks, 
at all other entries. 



used to evaluate k n at all other entries. 2N 
modes used. 5 panels were used to discretize the boundary. 



The FFT is 
1 is the total number of Fourier 



To test the accuracy of the approach, we have solved a both interior and exterior Dirichlet 
problems on each of the domains given in Figure [2] Exact solutions were generated by placing a 
few random point charges outside of the domain where the solution was calculated. The solution 
was evaluated at points defined on a sphere encompassing (or interior to) the boundary. The 
errors reported in Tables [3]j5] are relative errors measured in the /°°-norm, \\u € — u||oo/IM|ooj 
where u is the exact potential and u e is the potential obtained from the numerical solution. 
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For all geometries, 10 digits of accuracy has been obtained from a discretization involving 
a relatively small number of degrees of freedom, due to the rapid convergence of the Gaussian 
quadrature. This is especially advantageous, as the most expensive component of the algorithm 
is the construction of the linear systems, the majority of the cost being directly related to the 
number of panels used. Further, the number of Fourier modes required to obtain 10 digits 
of accuracy is on the order of 100 modes. Although not investigated here, the discretization 
technique naturally lends itself to nonuniform refinement in the rz-plane, allowing one to resolve 
features of the generating curve that require finer resolution. 

The number of correct digits obtained as the number of panels and number of Fourier modes 
increases eventually stalls. This is a result of a loss of precision in determining the kernels, as 
well as cancelation errors incurred when evaluating interactions between nearby points. This is 
especially prominent with the use of Gaussian quadratures, as points cluster near the ends of 
the panels. If more digits are required, high precision arithmetic can be employed in the setup 
phase of the algorithm. 



N P 


2iV+l 




25 


51 


101 


201 


401 


5 


3.9506E-04 


4.6172E-04 


4.6199E-04 


4.6203E-04 


4.6204E-04 


10 


1.3140E-05 


1.1091E-08 


4.8475E-09 


4.8480E-09 


4.8481E-09 


20 


1.7232E-05 


7.7964E-09 


4.7197E-12 


4.7232E-12 


4.7237E-12 


40 


2.7527E-05 


2.7147E-08 


2.8818E-14 


5.7173E-14 


5.7658E-14 


80 


2.118E-05 


9.4821E-09 


2.1529E-13 


2.0392E-13 


2.0356E-13 



Table 3. Error in internal Dirichlet problem solved on domain (a) in Figure^ 



N P 


2iV + l 




25 


51 


101 


201 


401 


5 


8.6992E-04 


1.3615E-03 


1.3620E-03 


1.3621E-03 


1.3621E-03 


10 


2.2610E-04 


9.6399E-05 


9.6751E-05 


9.6751E-05 


9.6752E-05 


20 


2.6291E-04 


4.6053E-07 


2.4794E-07 


2.4794E-07 


2.4794E-07 


40 


3.1714E-04 


2.8922E-07 


2.2875E-11 


2.3601E-11 


2.3605E-11 


80 


3.0404E-04 


3.5955E-07 


3.3708E-11 


3.3138E-11 


3.3150E-11 



Table 4. Error in external Dirichlet problem solved on domain (b) in Figure^ 



Np 


2iV + l 




25 


51 


101 


201 


401 


5 


4.3633E-04 


7.9169E-05 


7.8970E-05 


7.8970E-05 


7.8971E-05 


10 


3.9007E-04 


6.8504E-07 


2.0274E-08 


2.0272E-08 


2.0272E-08 


20 


3.8803E-04 


6.4014E-07 


3.2138E-11 


3.1624E-11 


3.1625E-11 


40 


3.8456E-04 


6.4098E-07 


6.5742E-12 


3.4529E-12 


3.4530E-12 


80 


3.9828E-04 


6.4486E-07 


6.8987E-12 


3.1914E-12 


3.1913E-12 



Table 5. Error in external Dirichlet problem solved on domain (c) in Figure [2j 

Finally, we investigated the conditioning of the numerical procedure. Figure [4] shows the 
smallest and largest singular values of the matrices {^1 + A( n )}^_ 200 on the domain shown in 
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1.5 






-200 



100 



200 



Figure 4. Maximum and minimum singular values for the matrices resulting 
from an 80 panel discretization of a sphere using 400 Fourier modes, where n is 
the the matrix associated with the n th Fourier mode. 



t igure |2j[a). The convergence of both the small est a nd the largest singular values to 1/2 follow 



from the convergence ||A( n )|| -4 0, cf. Section 2.2 Figure ^ indicates both that all matrices 



involved are well-conditioned, and that truncation of the Fourier series is generally safe. 
8.2. Helmholtz Equation. In this section, we repeat many of the experiments reported in 



Section 8.1, but now for the Helmholtz equation on an exterior domain with the associated 
"combined field" BIE formulation (461). The algorithm employed to solve the integral equation 
is the same as described in Section 4J with the caveat that the kernels are calculated using the 
fast procedure described in Section 



6.1 



Table [6] presents timing results, with all variables defined as in Section 8.1 



The largest problem size considered here has 80 panels and 201 Fourier modes, leading to 
160 800 unknowns discretizing the surface. Note that this problem size is slightly smaller than 



the largest considered in Section 8.1, due to memory constraints. This is because the matrices 
now contain complex entries, and thus use twice the memory compared with the Laplace case. 
The total running time of the algorithm is 97 seconds for this problem size. If given additional 
right hand sides, we can solve them in 0.51 seconds. We also remark that the asymptotic 
scaling of the cost of this algorithm is identical to that of the Laplace case, it simply takes more 
operations (roughly twice as many) to calculate the kernels and perform the required matrix 
operations. 

We have assessed the accuracy of the algorithm for various domains, discretization parameters, 
and Fourier modes. The boundary conditions used are determined by placing point charges 
inside the domains, and we evaluate the solution at random points placed on a sphere that 
encompasses the boundary. In the combined field BIE (46), we set the parameter v = k. 

First, we consider the ellipsoidal domain given in Figure [2^a). The major axis of this ellipse 
has a diameter of 2, and its minor axes have a diameter of 1/2. Table [7] lists the accuracy 
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N P 


2N + 1 


^mat 




T m 


^apply 


5 


25 


4.51E-02 


3.78E-03 


2.59E-04 


1.77E-04 


10 


25 


1.18E-01 


2.03E-02 


4.42E-04 


8.38E-04 


20 


25 


3.77E-01 


1.48E-01 


8.47E-04 


3.05E-03 


40 


25 


1.27E+00 


9.71E-01 


1.71E-03 


1.19E-02 


80 


25 


5.37E+00 


8.20E+00 


3.68E-03 


8.17E-02 


5 


51 


8.06E-02 


6.91E-03 


5.09E-04 


3.66E-04 


10 


51 


2.25E-01 


4.27E-02 


1.10E-03 


1.00E-03 


20 


51 


7.53E-01 


2.85E-01 


2.08E-03 


6.13E-03 


40 


51 


2.69E+00 


1.95E+00 


3.83E-03 


2.79E-02 


80 


51 


1.01E+01 


1.47E+01 


7.55E-03 


1.40E-01 


5 


101 


1.57E-01 


1.36E-02 


1.33E-03 


8.13E-04 


10 


101 


4.74E-01 


8.20E-02 


2.55E-03 


3.05E-03 


20 


101 


1.58E+00 


5.27E-01 


5.03E-03 


1.17E-02 


40 


101 


5.95E+00 


3.81E+00 


1.01E-02 


5.67E-02 


80 


101 


2.11E+01 


2.72E+01 


1.90E-02 


2.39E-01 


5 


201 


3.02E-01 


2.50E-02 


2.56E-03 


1.64E-03 


10 


201 


9.40E-01 


1.56E-01 


5.09E-03 


5.71E-03 


20 


201 


3.23E+00 


1.02E+00 


1.01E-02 


2.15E-02 


40 


201 


1.19E+01 


7.32E+00 


2.05E-02 


9.94E-02 


80 


201 


4.35E+01 


5.38E+01 


4.04E-02 


4.67E-01 


5 


401 


5.90E-01 


5.046E-02 


5.66E-03 


3.48E-03 


10 


401 


1.86E+00 


2.97E-01 


1.05E-02 


1.06E-02 


20 


401 


6.60E+00 


2.04E+00 


2.20E-02 


4.75E-02 


40 


401 


2.40E+01 


1.46E+01 


4.42E-02 


1.98E-01 



Table 6. Timing results in seconds performed for a spherical domain. 



N P 


2iV~+l 




25 


51 


101 


201 


401 


5 


4.2306E-04 


4.0187E-04 


4.0185E-04 


4.0185E-04 


4.0185E-04 


10 


3.4791E-06 


2.3738E-06 


2.3759E-06 


2.3760E-06 


2.376E-06 


20 


7.8645E-06 


4.5707E-09 


7.1732E-11 


7.1730E-11 


7.173E-11 


40 


1.3908E-05 


1.5980E-08 


2.9812E-13 


3.1500E-13 


3.148E-13 


80 


1.0164E-05 


5.4190E-09 


4.9276E-13 


4.8895E-13 





Table 7. Relative error in external Helmholtz problem for the domain in Figure 
p[a). The domain is 1 wavelength in length (the major axis). 



achieved. We achieve 9 digits of accuracy in this problem with 20 panels and 51 Fourier modes. 
We have not padded the two sequences in the convolution procedure described in Section 6T but 
as more modes are used the tails of the sequence rapidly approach zero, increasing the accuracy 
of the convolution algorithm utilized to calculate the kernels. Table [8] displays the same data, 
but with the wavenumber increased so that there are 25 wavelengths along the length of the 
ellipsoid. We see a minor decrease in accuracy as we would expect, but it only takes 40 panels 
and 51 Fourier modes to achieve 8 digits of accuracy. 

We now consider the more complex domains given in Figures ^h) and[2](c). They are a wavy 
shaped block and a starfish shaped block with the outer diameter of size roughly 1.5 and 1.0, 
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iV P 


2iV + l 




25 


51 


101 


201 


401 


5 


2.0516E+00 


3.4665E+00 


4.1762E+00 


4.4320E+00 


4.4951E+00 


10 


2.3847E-01 


2.4301E-01 


2.4310E-01 


2.4312E-01 


2.4313E-01 


20 


1.7792E-02 


8.8147E-06 


8.8075E-06 


8.8081E-06 


8.8082E-06 


40 


1.7054E-02 


5.4967E-08 


3.7994E-10 


3.7999E-10 


3.7998E-10 


80 


1.7302E-02 


1.6689E-08 


1.8777E-11 


1.8782E-11 





Table 8. Relative error in external Helmholtz problem for the domain in Figure 
pi a). The domain is 25 wavelengths in length (the major axis). 



N P 


2N + 1 




25 


51 


101 


201 


401 


5 


2.6874E+00 


2.5719E+00 


2.5826E+00 


2.2374E+00 


1.9047E+00 


10 


1.6351E+00 


4.2972E+00 


4.6017E+00 


1.0380E+00 


1.0328E+00 


20 


4.3666E+00 


2.6963E-03 


2.6966E-03 


2.6967E-03 


2.6967E-03 


40 


4.4498E+00 


9.1729E-08 


4.1650E-08 


4.1661E-08 


4.1664E-08 


80 


4.3692E+00 


7.3799E-08 


1.4811E-10 


1.4812E-10 





Table 9. Relative error in external Helmholtz problem for the domain in Figure 
pfb). The domain is 10 wavelengths in length (the major axis). 



N P 


2iV + l 




25 


51 


101 


201 


401 


5 


1.1140E+01 


3.2400E+01 


3.3885E+01 


4.4322E+01 


1.6151E+01 


10 


1.9626E+01 


7.7739E+01 


6.2931E+01 


3.6335E-01 


3.6344E-01 


20 


2.6155E+01 


4.9485E+01 


2.5796E+01 


4.8239E-05 


4.8245E-05 


40 


3.3316E+01 


5.0645E+01 


2.4330E+01 


1.3841E-09 


1.3846E-09 


80 


1.6966E+01 


6.0163E+01 


2.4354E+01 


1.6510E-10 





Table 10. Relative error in external Helmholtz problem for the domain in Figure 
pfc). The domain is 10 wavelengths in length (the major axis). 



respectively. The accuracy for various values of Np and N are given in Table [9] and 10 We 



achieve 8 digits of accuracy with 40 panels and 51 Fourier modes for the wavy block and 9 digits 
of accuracy with 40 panels and 201 Fourier modes for the starfish block. 
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